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Abstract 

The use of mathematical models has helped to shed light on countless phenomena 
in chemistry and biology. Often, though, one finds that systems of interest in these 
fields are dauntingly complex. In this paper, we attempt to synthesize and expand 
upon the body of mathematical results pertaining to the theory of multiple equilibria 
in chemical reaction networks (CRNs), which has yielded surprising insights with 
minimal computational effort. Our central focus is a recent, cycle-based theorem by 
Gheorghe Craciun and Martin Feinberg, which is of significant interest in its own 
right and also serves, in a somewhat restated form, as the basis for a number of 
fruitful connections among related results. 

1 Introduction 

Chemical reactions have long been a popular and fruitful subject for exploration with 
mathematical models. A wide range of behaviors can emerge from the evolution of a 
chemical system over time, and these can frequently be explained in a useful way by 
analyzing a dynamical model. Even in complex biological systems, one often finds it 
possible to isolate tractable networks of interacting proteins and to shed light on certain 
of their fundamental properties using fairly simple chemistry and mathematics. 

In this paper, we will focus on the capacity or incapacity of chemical reaction networks 
(CRNs) to admit multiple equilibria in their dynamics. Naturally, the topic of multista- 
bility (or its opposite, injectivity, depending on one's perspective) has attracted interest 
over the years in the chemistry, physical chemistry, and chemical engineering literature 
[H H21 H31 [151 [161 [HI EU [22I ^ ST]. More recently, biologists have also come to recog- 
nize the importance of multistability for living systems, especially on the molecular scale, 
where the presence of multiple equilibria can be associated with switch-like processes and 
biological memory, particularly in cell differentiation and development. Specific areas 
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of study have included gene networks [10J [251 ED EQ], enzyme catalysis and induction 
PH Q2J OH [39], and cellular signaling [H [231 [2ZI [281 E31 EH1 SS]. A number of the re- 
sults we will discuss have also come from farther afield, with applications ranging from 
ecology [2] to economics [HI EH EH HE] to probability theory [33] . 

In describing a CRN as a dynamical system, we will be dealing with the ways in which 
some state variables - the concentrations of chemical species in the system - influence the 
rates of change of others over time. One species can activate or inhibit a second species, 
which can in turn have an effect on a third, and so on. In the study of multistability, we 
will be particularly concerned with the (usually indirect) influence that one state variable 
exerts on its own value. This kind of "cycle" (as we will call it) is often easy to describe 
qualitatively, simply by examining the reactions in the network, and in many cases this 
proves to be enough to yield important and subtle information about the properties of 
the network [IB] . For example, if we were to observe multistability in a chemical process 
composed of several unknown steps, we could apply the results in this paper to help 
us decide among proposed alternatives for the underlying series of steps, even if the 
possibilities were superficially very closely related. 

A recent theorem by Gheorghe Craciun and Martin Feinberg [151 CSj is the central 
focus of our paper. Over the last twenty years, Feinberg has been the leading developer 
of "chemical reaction network theory," or CRNT. His earlier work on multistability was 
based on a property known as the deficiency of a network [211 E21 [35], but more recently, 
he and his colleagues have focused on cycle-based conditions [T5l [161 HB1 SZ] , which can 
be compared and combined with other similar results in a number of productive waysQ 
Feinberg has not been widely cited in either the biology or the mathematics literature, 
so one of our goals is to bring attention to what has been an underappreciated body of 
work. 

We begin with a series of definitions and a discussion of some of their implications in 
sections [2] - [31 followed by some general results in section H] regarding global injectivity. 
Section [5] introduces graphs and cycles and the Thomas-Soule theorem, and section [6] 
presents a generalized and substantially reformulated exposition of Craciun-Feinberg. We 
conclude with a discussion of the relationships among the various Jacobian-based theorems 
on multistability in section [7] and an example of the theory in action in section [HI 

2 Preliminary definitions and notation 

Imagine a container, either a laboratory vessel or a living cell, within which some chem- 
ical reactions take place at constant temperature. Let the active chemical species be 
Si, S 2 , ■ ■ ■ , S n , and suppose there are m reactions that occur among them. The j th reac- 
tion is given by CjiSi + c j2 S 2 + h c jn S n —> djiSi + d j2 S 2 + h d jn S n , where the 

and dji are nonnegative integers. Some reactions may be reversible, in which case we use 

1 While deficiency theory is certainly a worthy subject for study, the mathematics behind it is fairly 
specialized, and we will not be dealing with the subject here. Similarly, we will not delve into the 
stoichiometric network analysis (SNA) of Bruce Clarke [TJ [T^l H3] ■ 
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the shorthand Cj\S\ + Cj 2 S 2 H — • + Cj n S n ±=t dj\S\ + dj 2 S 2 + • — h dj n S n . For now, we will 
consider the forward and reverse directions to be two separate reactions. 

At time t, we write Xi(t) for the molar concentration of species Si in the container, 
with Xi(t) G M + for each i and all t, and let x be the vector (xi, x 2 , ■ ■ ■ , x n ). (At times, we 
will also use the notation [Si] to denote the concentration of Si. These quantities change 
over time due to the reactions that take place, which we will assume are elementary and 
are thus governed by mass-action kinetics (see section [372]) . This means that the rate (in 
concentration per time) of reaction j at a certain point in time is given 
where kj is a positive constant, referred to as the rate constant for the reaction. From 
stoichiometry, then, we know that the rate of change of Xi due to the effects of reaction j 
is (dji — Cjijkj niLi x< i Jl ■ T° simplify the notation, let us write Cj = (cji, Cj 2 ■ ■ ■ , Cj n ), 
dj = (dji, dj 2 . . . , dj n ), and u v = n a M a Q f° r arbitrary vectors u = (u a ) and v = (v a ). 
Summing over all reactions yields the complete rate equation 

d m 

—jj- = ^^(dji — Cjijkjyi J . 

3=1 

Note that a reaction will only contribute to this sum when and dji are different. As 
a special case, if 5^ does not participate in reaction j, then Cji = dji = 0, and reaction j 
does not affect the concentration of Si, as we should expect. Henceforth, we will make 
the assumption that in fact Cji and dji cannot both be positive, i.e. no species appears on 
both sides of the same reaction. To make the notation cleaner, we will write = dji — Cji 
and ej = (e^i, tj 2 . . . , e jn ), remembering that either dji or Cji must be zero. In terms of 
chemistry, this means that we require the processes inside the container to be broken down 
into sufficiently simple mechanisms so that each reactant in a given reaction undergoes a 
chemical changell 

The changes in the chemical species over time can be condensed into a single function 
as follows: 

Definition 1. The vector rate function for a chemical system is the polynomial func- 
tion F(x) from (M+) n to W 1 defined by F(x) = 

From F, we define the associated Jacobian matrix Jp = (Jik), where 

j _ <9Fi _ _d_ ( dxj_ \ _ c^ e fc ^ 
lk dx k dx k \dt 4-1 x k 31 3 

V / j=l 

We shall see in section 13.41 that it is valid to assume that the Xi remain positive for all t, 
so it is not a problem for them to appear in the denominators. Denote by Jijk the j th 
term in Jik- 

Jijk CjikjU. ^ . 

This quantity represents the effect of on the concentration of Si via reaction j. 
2 For further discussion, including the issue of catalysis, see section [531 
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3 Remarks on the definitions 



The framework introduced in section [2] contains several assumptions with important math- 
ematical implications. Here we will briefly discuss the chemical and mathematical ratio- 
nale for several of our premises. 

3.1 Uniformity and continuity 

We will assume that our reaction vessel is isothermal and spatially homogeneous, and 
that the reactions take place in the liquid phase at constant density. These conditions 
ensure that the only time-dependent quantities in the system are the concentrations of 
the reacting species. In practice, depending on what type of container we are working 
with, the conditions may be not be perfectly satisfied: for example, a living cell may be 
close to isothermal, but it is far from well-mixed, although we often do not know enough 
to justify more precise assumptions. Our hope is that networks of interest take place 
within a sufficiently uniform environment so that the qualitative behaviors of interest are 
unaffected by disregarding the heterogeneities. 

A related assumption we will make is that the Xi are continuous (in fact, differen- 
tiable) functions of time. In real life, molecules exist in discrete units, and cells may have 
extremely small numbers of certain proteins, meaning that stochastic modeling can some- 
times be more appropriate [ST] . However, it is standard to assume that concentrations are 
continuous - especially in the lab - and this will be an essential condition for the results we 
present. Once we assume continuity, differentiability follows naturally from the explicit 
rate functions. 

3.2 Mass-action kinetics 

As stated in section [2j we will be assuming that all of our reactions are written at the 
elementary level, where the law of mass action kinetics makes good physical sense. Take, 
as an example, an elementary reaction A + B — > C '. One molecule of A will combine with 
one B to form a C whenever the two reactant molecules collide with a certain energy 
and orientation. If we assume that the probability of a collision is proportional to the 
abundances of the two species, then we immediately have the expression fc[A][J3] for the 
reaction rate, where k is a constant. This is only an approximation, and in certain physical 
situations it will be more valid than others, but it is almost universally accepted |34j . 

Most of the time, especially in biology, an observer will see a composite process taking 
place rather than the mechanisms underlying it (see section [3T31 for a basic example), and 
it can be very difficult to write down all of the elementary reactions involved in a CRN. 
The results in this paper will often be useful in discriminating among several proposed 
mechanisms for a process, based on the behaviors that each would permit according to 
the theory [4~7] . 

Mathematically speaking, the most important consequence of the law of mass action 
is that vector rate functions arising from CRNs will always be polynomials. 
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3.3 Catalysis and the positions of species in reactions 

We assumed in section |2] that no species appears on both sides of the same reaction. Our 
claim here is that any proposed reaction with a species as both a reactant and a product 
can be rearranged or decomposed into a simpler series of steps. 

As an example, suppose we had a proposed reaction A + B — > 2A + 2C. It could be 
that one molecule of A is a spectator and does not participate chemically; in this case, 
the true reaction mechanism is B — > A + 2C. It could also be that the molecule of A 
on the reactant side undergoes a preliminary, implicit chemical change, so that the true 
mechanism consists of multiple steps, perhaps A + B — > D + 2C followed by D — > 2A 

We note in particular that the conversion of a substrate S into a product P through 
the action of an enzyme or other catalyst E can be decomposed into multiple reactions 
in this way. The total effect of such reactions is E + S — ► E + P, but the rates that 
are observed in catalysis are not of the form /cfPjf.S'], so there must be some intermediate 
steps involved. The most common model is the Michaelis-Menten mechanism, which 
involves a reversible enzyme-substrate-complex step: E + S ES — > E + pj§ Assuming 
that the concentration of E is very small compared to that of S, the full reaction rate 
(i.e. for the formation of P) under this model is given by ki[S]/(k2 + [S]), where h\ and 
hi are constants, and this expression matches observations quite well [H]- Thus, while 
the process begins and ends with one molecule of E present, in none of the elementary 
reactions does E or any other species appear on both sides. 

Some of the best-known instances of multistability in biochemistry involve enzyme- 
catalyzed reactions, particularly when several enzymatic modules are coupled [16J, I2Z] , as 
we will see, for example, in section [8j 



3.4 The species domain 

It is clear that in the systems we are considering, the state variables, which represent 
concentrations of chemical species, cannot be negative. In fact, the assumptions of conti- 
nuity and mass-action kinetics allow us to restrict the domain to (M + ) n , as we stated in 
section El this will be useful to us at several points in the discussion|3 

If Xi is initially zero for some % (say at time t — 0) and Si is not produced via any of the 
reactions in our network, then Xi will always be zero, and we can simply ignore Si and all 
reactions in which it partcipates. Otherwise, assume that x.;(0) > and that Xi(to) = 
for some to > 0, where we are given some initial conditions such that the Xj have finite 
solutions through time t . We can choose t such that Xi(t — e) > for any e > 0. As in 
section [21 we have dxi/dt = YlT=i e jikjX° s - Since Si is assumed never to appear on both 
sides of the same reaction, we know that < if and only if x c J contains a factor of Xi, 
as these two conditions hold exactly when Si is on the reactant side of reaction j. So, 
we can write dxi/dt = — Xi ■ Pi(t) + P2{t), where Pi and P2 are polynomials with positive 

3 In fact, the reaction ES — > E + P, like most chemical reactions, will probably be slightly reversible, 
but given the energetics of catalysis, it is almost always assumed to be unidirectional, especially when 
the product P is being removed from the area of formation [51] . 

4 Craciun and Fcinbcrg make this assumption as well, although without including a detailed analysis. 
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coefficients in terms of the x\ (and thus implicit nonnegative functions of t) and P2 does 
not contain any factors of If < e < t' < to, then dividing by Xi and integrating 
yields, up to a constant, 




Pi(t) + 




dt. 



As t' — > to, we have Xi — > + , and so lnxj — > —00. This means that Pi(t') — * 00, in 
order to make the right-hand side approach —00, and hence the concentration of some 
species goes to infinity in finite time. This is a contradiction, so we conclude that the 
concentration of must remain strictly positive for as long as the system has solutions. 

4 Derivatives and conditions for injectivity 

In chemistry or thermodynamics, or even in everyday usage, an equilibium is a state 
of balance among the constituents of a system. We will be interested in studying the 
equilibria (also referred to as the steady states) of a reaction network, specifically to 
determine whether or not there may be more than one. To be precise, an equilibrium 
point is a concentration vector x such that F(x) = oj^l Equivalently, if a CRN starts at 
an equilibrium point and evolves in time, it will not move from this point, as any change 
in a concentration due to some reaction is balanced by the other reactions in the network. 

We will say that a network and its rate function are injective if F(xi) 7^ F(x 2 ) 
whenever x x 7^ x 2 . This is the standard definition of an injective function. Note that an 
injective CRN cannot have mulitiple equilibria. 

For a two-dimensional dynamical system dx/dt = f(x,y) and dy/dt = g(x,y), the 
easiest way to check for equilibria is usually to graph the nullclines f(x,y) = and 
g{x,y) = and to find their intersection points. The stability properties of these steady 
states will also be fairly easy to determine^ In real life, however, and especially in 
biochemistry, most CRNs of interest will have more than a few species. With polynomial 
rate functions, there might be some hope of using tools from algebraic geometry, but in 
any case, the set of equations will be intimidating. In the next few sections of this paper, 
we will begin to consider the question of how to determine, with minimal effort, when a 
given network is injective, with derivatives as our primary tool. 

4.1 Simple examples 

In one dimension, if / : R — > R is a differentiable function and f(a) = f(b) with a < b, then 
Rolle's theorem (or the mean value theorem) tells us that /'(c) = for some a < c < b. 
Thus, if f'(x) is nowhere zero, then f(x) is injective. This is about as simple as the 

5 The notion of steady states is sometimes taken to include limit cycles, but we will restrict the 
definition here to points. 

6 Stability is a very rich and important problem and one that is obviously relevant to biology, but as 
such it deserves much more space than we could devote to it here. 
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relationship between injectivity and zero-valued derivatives can be; in higher dimensions, 
the story becomes much more complicated. 

Suppose F is a differentiable, vector-valued function from IR n to W 1 . The "derivative" 
of F is now a matrix, the Jacobian = (Jik), where Jn~ = dFi/dxk- By the inverse 
function theorem, F is locally injective in the neighborhood of points where det(J^) is 
nonzero, but what can we say about its global properties? In [23], Gale and Nikaido pro- 
vide the following example: let F : R 2 — > 1R 2 be defined by F(x, y) = (G(x, y), H(x, y)), 
where G(x, y) = e 2x — y 2 + 3 and H(x, y) = 4e 2x y — y 3 . Then, 



However, F(0, 2) = F(0,—2) = (0,0). Thus, even in two dimensions, an everywhere 
nonzero Jacobian determinant does not imply global injectivity. In fact, this function was 
given in 1965 as a counterexample to a weaker conjecture by Samuelson [36], which stated 
that F is injective if all of the leading principal minors of - i.e., the determinants of 
the upper-left ixi submatrices - are nonvanishing. So, any theorem of this sort will need 
to be weaker still. 

4.2 The Gale-Nikaido theorem 

Having found a counterexample to Samuelson's conjecture, Gale and Nikaido were able 
to formulate and prove two stronger conditions for global injectivity based on the Jaco- 
bian Jp [23] • Their work has subsequently inspired a number of extensions [43J, but the 
original Gale-Nikaido theorem is all that will interest us here. 

We will say that M is a P-matrix if all of its principal minors (the determinants of the 
submatrices consisting of all entries with both indices in some given subset of {1, 2, . . . , n}) 
are strictly positive and that it is a weak P-matrix if the determinant of M is positive 
and all other principal minors are nonnegative. 

Theorem 1 (Gale-Nikaido [24]). Let D be a rectangular region ofR n , and let F : D — > M n 

be a differentiable function. If J p is a P-matrix at all points x £ D, then F is injective. 
The same result holds if Jp is a weak P-matrix for all x £ D, provided that D is open. 

The proof for P-matrices is based on linear inequalities. In particular, if we write 
y ^ z to mean that ?/j > for each component i, it is shown first that the only solution 
to the system J^y ^ 0, y ^ is y = 0. Then, more generally, for any z £ D, the system 
F(y) ^ -^(z), y ^ z is only satisfied (in D) by y = z. 

The proof for weak P-matrices, meanwhile, involves several ingredients from topolog- 
ical degree theory (also known as index theory). First, if z is some point in D, then the 
function G(y) = F(z) + y — z clearly only takes the value F(z) at y = z. Also, the 
Jacobian of the function F(y) = F(y) + Ay for A > is a P-matrix, so F is injective. 
This allows us to construct a homotopy from F to G via F, and hence the degrees of 



=>■ det(Jjr) 
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boundary cycles of D are the same under F and G. The last useful fact is that these 
degrees equal the sums of the degrees of all points y G D with F(y) = F(z) (respectively 
for C) 13 Since det(Jp) is positive and det(Jc) clearly is as well, all of these degrees are 1. 
Thus, because G is injective on D, so is F. 

4.3 Polynomial maps and the real Jacob ian conjecture 

The results we have seen so far are very general, applying to all differentiable functions. 
The vector rate function for a CRN, however, has a specific form. First and foremost, if 
we assume mass-action kinetics, then each component of the rate function is a polynomial 
in terms of the concentrations of the n chemical species. Because polynomials tend to 
be particularly well behaved, we might hope for a more powerful injectivity theorem 
in this special casej^l Indeed, a great deal of effort has been devoted to the "Jacobian 
conjecture" (for complex functions) and the related "real Jacobian conjecture," which 
posits that a polynomial function F : IR n — > IR n with nonvanishing Jacobian determinant 
is globally injective. In 1994, however, Pinchuk proved that we are not so lucky, providing 
a counterexample, a polynomial function F = (p, q) : IR 2 — > M 2 with deg(p) = 10 and 
deg(g) = 35 that has a nowhere zero Jacobian determinant but is not injective |41j . 

With the real Jacobian conjecture shown to be false, others have attempted to prove 
a modified version (see [20J for a thorough reference). One easy case is when all of the 
components of F are first-degree polynomials: then F is a linear transformation, and 
it is well known that F is injective if and only if is nonsingular. If F is quadratic, 
the conjecture still holds [53], and the proof is quite simple [ID]. Beyond this, however, 
the theory becomes much harder; an example is the following, by Jelonek [30]. For a 
continuous map F : X — > Y, we say that F is proper at a point y G Y if there exists 
a neighborhood U of y such that F^iU) is compact. Let Sf be the set of points in Y 
at which F is not proper. Jelonek's theorem states that if F : W l — > W 1 is a polynomial 
map with nowhere zero Jacobian determinant and Sp has codimension at least 3, then F 
is injective. 

This theorem and others like it, while providing a variety of potentially useful results 
on injectivity, tend to include conditions that are very difficult to apply to CRNs, for 
example in relation to the properties of the set Sf- It is also significant that they are 
based on the requirement that det(J.p) 7^ for all x G M n (or C n ), whereas our domain 
of interest is only (BL + ) n . 

5 Graphical conditions and Thomas-Soule 

So far, we have seen the important role that the Jacobian of the vector rate function 
plays in determining the number of equilibria of a CRN. In practice, though, the Jacobian 

7 As observed in [51], this theorem, which is cited as the "Kronecker theorem on indices" in [24), is 
remarkably similar to Gauss's law in electrostatics. 

8 For example, Campbell [9j has shown that Samuelson's conjecture holds for polynomial functions 
F : R n -> R n . 
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determinant will be an immensely complicated polynomial, and finding its zeros will be 
a daunting task. Might there be a way to use some simpler properties of the chemical 
system to determine, with limited effort, whether or not the Jacobian determinant can 
vanish in a certain domain and with certain rate constants? Naturally, we expect any 
such conditions to sacrifice some generality for the sake of elegance. In the sections that 
follow, we examine some approaches based on a network's interaction graph. 

5.1 The interaction graph and its cycles 

Definition 2. The interaction graph G(x) of a CRN is the oriented graph with one vertex 
for each species Si and an edge from S^ to S i2 {i\ ^ i 2 ) if and only if J i2il is nonzero. 
Each edge is given a sign, which is the sign of J i2il . 

Note that the signs of the edges in the interaction graph can change for different con- 
centrations x. We will refer to an edge for which this occurs as concentration-ambiguous. 
A concentration- ambiguous edge from S^ to S i2 will disappear when J iail = 0, but be- 
cause our results are stated in terms of the possible configurations of G(x) for any x, this 
will not matter. 

Although the diagonal terms Ju of the Jacobian are usually nonzero, we omit self-edges 
from the graph because they are not used to form cycles (Definition [3]). This is another 
way of saying that autocatalysis and autoinhibition are expressed indirectly at the mass- 
action level (as in section 13. 31) . As we will see, diagonal terms are always nonpositive 
(and, under certain assuptions, always negative). 

Definition 3. A cycle in an interaction graph is an ordered subset (z 1; i 2 , . . . , in) of 
(1, 2, . . . , n) with 1 < N < n such that there exists an edge from S^ to Si 2 , from Si a to 
Si 3 , and so on, including fron Si N to S^. A cycle is classified as positive or negative based 
on the product of the signs of all of its edges. 

In some examples, when the species names are more prominent than the indices, we 
will also use the notation Si 2 =$>■■•• => S.^ to indicate a cycle. Two cycles are 

disjoint if they have no vertices in common; otherwise, they intersect. 

We will be interested not only in the sign of a cycle but also in its sub-sign. For an 
edge from S^ to S i2 , we define its sub-sign, for some chosen reaction j that contributes 
a term Ji 2 j^ to Ji 2 i ± , to be the sign of J% 2 ji x . The sub-sign of a cycle is then the product 
of the sub-signs of its edges, given N chosen reactions: sub-positive if there are an even 
number of sub-negative edges and sub-negative if there are an odd number of sub-negative 
edges. The sub-sign of an edge could depend on the choice of j, so we will always need 
to specify which reaction indices we are using. If the edge is sub-positive for some j and 
sub- negative for some other j, we refer to it as reaction- ambiguous. An edge could be both 
concentration-ambiguous and reaction-ambiguous (see in particular Lemma EJ). However, 
because all of the concentrations Xi are positive, the sign of Ji 2 ji 1 does not depend on x. 

Given two cycles and choices of reactions j as above for each of their edges, we will say 
that the cycles strongly intersect if they have at least one vertex Si in common such that 
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the reaction j that is chosen for the outgoing edge from Si is the same for both cycles. If 
this condition is not met, then they are weakly disjoint. 

5.2 An example 

To illustrate our definitions, consider the following hypothetical network, to which we will 
return in section 17721 

1. A + B + C -> X 

2. A + B + D 

3. C + E->A 

4. D + E -> B 

5. A -> Z 

6. Z^D. 

Using the law of mass action, we can compute the rate functions for each species: 

-h[A][B][C] - k 2 [A] [B] [D] + h[C][E] - k 5 [A] 
-h[A][B][C] - k 2 [A][B][D] + h[D][E] 
-h[A][B][C] - h[C][E] 
-k 2 [A][B][D]-k 4 [D][E] + k 6 [Z] 
-k 3 [C][E]-k 4 [D][E] 
h[A][B][C] 
k 2 [A][B][D] 
h[A]-k G [Z]. 

The interaction graph for this network contains many cycles; one of interest is C A 
Z =^> D B =^> C, which we will call cycle Numbering the species in alphebetical 
order from 1 to 8, the entries in the Jacobian which comprise this cycle are J13, J$i, J48, J24, 
and J32, which have the values 

9 An example of a cycle intersecting c is d : D B ^> D. If we pick the term J224 in c and the terms 
J244 and J422 in d , then this is not a strong intersection. 



d[A\ 

dt 
d[B] 

dt 

d[C] 

dt 
d\D] 

dt 

dt 
d[X) 
dt 

d\Y] 
dt 

dt 
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h 
h 

-k 2 [A][B} + h[E] 
-h[A][C}. 

Let us suppose for later use that k\ = k 2 and k% = k±, remembering also that all 
of the kj are positive. Using the terms Ji 33 , J 85 i, J 468 , J224, and J 3 i 2 , c is sub-positive; 
the first three terms (A; 3 [i?], k 5 , and k s ) are positive, while the last two terms (— &2[A][-B] 
and — ki[A] [C]) are negative. However, J 13 and J24 are identical, so their signs are always 
equal, and hence their product is nonnegative and c is always negative. 

5.3 The Thomas- Soule theorem 

In 1981, Rene Thomas conjectured that any dynamical system displaying stable oscil- 
lations must have at least one negative cycle, while any system with multiple steady 
states must contain a positive cycle [52J. Intuitively, these associations of behaviors 
with signed cycles make sense. For example, suppose we have a three-species cycle 
A B =>- C =>- A. If all three edges are negative, then the cycle is negative. 

In this case, if the concentration of A is increasing, it will cause the concentration of B to 
decrease, which will in turn cause the concentration of C to increase, which will cause the 
concentration of A to decrease. Since we initially assumed A to be trending upward, we 
can see how the cycle promotes oscillatory behavior. By contrast, if one edge is positive 
and the other two are negative, then the cycle is positive, and an increase in A is self- 
reinforcing. So, if the concentration of A is perturbed away from equilibrium, the positive 
cycle will aid in pushing the system into a new basin of attraction. 

In the years after Thomas's paper, a number of authors found the conjecture to be true 
in special cases [HI [261 H21 HE], and in 2003, Christophe Soule presented a full proof [19]. 
While the Thomas-Soule theorem is certainly elegant, we should note that it gives no 
information about the sufficiency of positive cycles: in order for a system to display 
multistability, certain other conditions will need to be met, either in the structure of the 
network or in the ranges of values of specific state variables or parameters. 

The precise statement of Soule's result is as follows. 

Theorem 2 (Soule [49]). Let D be an open, rectangular region ofW 1 , and let F : D — > R n 

be a differentiate function. If the interaction graph of F has no positive cycles for any 
x e D, then F is injective. 
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The general idea behind the proof is that each term in a given principal minor with TV 
elements corresponds to a permutation a of {1, 2, . . . , N} such that the (i, a(i)) th entry 
in the submatrix is nonzero. This a in turn corresponds, by taking its algebraic cycle 
decomposition, to a union of disjoint cycles in our sense of the word (some possibly of 
length 1) that together incorporate all N elements. As we will see in section 16.21 if all 
of the cycles in the interaction graph are negative, then all of the terms in the expansion 
of any principal minor will have the same sign, namely (— 1)^. Neither the link between 
cycles and determinants [H] nor this result about negative cycles [6] are new in Soule's 
work, but these facts are typically used in connection with the full Jacobian. Because 
nonsingularity does not by itself imply injectivity, however, it is important that Soule 
takes the results a step farther by involving the principal minors of J p. 

If we apply the above reasoning about cycles and determinants to the matrix — Jp, we 
see that if G(x) never contains a positive cycle, then none of the minors of this matrix 
will ever be negative. If the full Jacobian determinant is strictly positive, then — Jp is 
a weak P- matrix, and by Gale-Nikaido (Theorem [T]), — F (and hence F) is injective. If 
we assume homogeneous outflow (see section [6]), then all minors have a nonzero diagonal 
term, and — Jp is in fact a P-matrix, meaning that the domain D need not be open. If 
neither of these assumptions is met, then the theorem still holds, with the caveat that it 
only guarantees a single nondegenerate zero, i.e. one at which det(Jp) ^ [49, 50J. 

Note that this theorem, like that of Gale-Nikaido, applies to any differentiable function. 
As a special case, of course, it can be applied to vector rate functions for reaction networks. 
However, some of its power is lost when we confine our attention to this restricted class 
of polynomials. The next theorem we will examine, by contrast, is formulated specifically 
for reaction networks. 

6 The Craciun-Feinberg theorem 

Recently, writing from the perspective of chemical engineering, Gheorghe Craciun and 
Martin Feinberg have presented an exciting, graph-based theorem on multistability in 
CRNs [H]-[I7]. They imagine their reactions to be taking place in what is called a 
continuous-flow stirred tank reactor (CFSTR), an isothermal, spatially uniform container 
with constant liquid flow streams in and out. If r is the rate of flow in units of 1/time, 
then the inward stream, or feed, is a constant vector rf (written with one component per 
species), with f in units of moles per volume, while the outward stream, or outflow, is 
assumed to be drawn homogeneously from the reactor and hence equals rx. We can think 
of the feed and the outflow as reactions of the form — > S and S — > 0, respectively, to 
be appended to the CRN, with the outflow, we note, obeying mass action kinetics. For 
example, in the network from section 15. 2\ the rate function for A would become 

M = rh-kMWm-HAmm + HCm-k^-riA}. 
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A network with this type of outflow is typically referred to as open or homogeneous. 

Upon computing any partial derivatives of the rate expressions, the feed term disap- 
pears, so the feed does not affect the Jacobian. From now on, we will ignore the feed 
stream, noting that because it is constant, the feed-free reaction network is injective if 
and only if the full CFSTR network is injective [15]. On the other hand, the outflow 
contributes to the diagonal terms of the Jacobian matrix, and this fact is used in the 
proof. One might argue that the requirement of outflow reactions makes the theory less 
widely applicable, but a few points may be offered in response. First, in a biological 
context, the outflow may be interpreted as the degradation of proteins or RNAs, which, 
as with the flow from a CFSTR, is typically assumed to proceed at a rate proportional 
to the concentrations of the species in question [101 [331 [50] Additionally, even if we 
wish to insist that the concentrations of certain enzymes or other species are not subject 
to outflow or degradation, Craciun and Feinberg prove that their results carry over, with 
only degenerate exceptions, to this "entrapped species" case [T7] . 

In addition to the presence of outflow, the Craciun-Feinberg theory differs from other 
approaches in the way that it treats rate constants and in how it parses the Jacobian 
matrix. So far, all of our results have been formulated to apply to a single, well-defined 
network, complete with rate constants. In Craciun-Feinberg, however, when one is able 
to prove injectivity for a CRN, it applies for any positive values of the kj. This added 
power can be useful, especially in biology, because accurate values for rate constants 
and concentrations can be very difficult to obtain, and realistic reaction networks will 
usually contain many such parameters [5j EES] . On the other hand, though, we might 
want to be able to distinguish, for a given network, which parameter regimes can support 
multistability and which cannot. The structure of the proof in section 16.31 is very much 
related to the parameter-independent nature of the theory. 

With regard to the Jacobian, while Craciun and Feinberg use different terminology 
from ours, their key innovation is effectively to split each entry into its component 
terms Jijk and to expand det(Jp) in terms of these components. This allows for finer 
criteria in relating the structure of G(x) to the value of det(Jp); some (sub-) positive 
cycles can be accounted for without violating the conclusion of injectivity. Section 16.21 is 
devoted to describing how this process works. 

6.1 New notation and the statment of the theorem 

Craciun and Feinberg's statements and proofs are somewhat different from ours, both in 
choice of language and in some more significant structural ways [15*1 [To] . As a result, we 
will begin by explaining a few of their concepts, in order to illuminate the relationships 
between their work and the more standard interaction-graph-based literature (about which 
we will say more in section [7]) . 

First, Craciun and Feinberg define two new graphs, called the SR and OSR Graphs, 
which have vertices representing both reactions and species as well as more elaborate edge 

10 As Craciun and Feinberg observe, their results are unchanged when the rate constants for the outflow 
are different for the different species |15j . 
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labels, and they make use of the cycles in these graphs (SR in the statement of the theorem 
and OSR in the proof) rather than those in the standard interaction graph. This discrep- 
ancy is the reason for our introduction of sub-signs in the interaction graph. Second, they 
classify cycles on the basis of parity, rather than sign, as follows. Let (ii, i 2 , . . . , iff) be a 
cycle in G(x), and let j\,j 2 , • • • , j/v be reaction indices such that Ji h+1 j h i h is nonzero for 
1 < h < N (cyclically). Each index jh is called shared if Si h and Si h+1 are on the same 
side of reaction j^. Then, we say that i 2 , . . . , ijv) has even potential if, for some choice 
of the indices ji,j 2 , ■ ■ ■ 1JN1 an even number of them are shared. For example, in our use 
of the cycle c in section 15.21 the shared edges are those from D to B and from B to C, 
so c has even potential. 

Having given a sense of their framework, we will now show how their definitions relate 
to the notions we have already introduced. Let {ii,i 2 , ■ ■ ■ , i/v) an d ji,j2, ■ ■ ■ , Jjv be as 
above, so that, for each h, 

j. . . r J" ; " , . . /.. x <-.i,, 

is nonzero. This implies that Cj h i h is nonzero, which is to say that Si h is on the reactant 
side of reaction j h . The sign of Ji h+1 j h i h is determined by that of ej hih+1 , so this term is 
positive when S ih+1 is on the product side of reaction j h and negative when S ih+1 is on the 
reactant side. (Recall that a species cannot appear on both sides of a reaction.) Thus, 
the edge Ji h+1 i h is sub-negative if and only if jh is shared. By counting the number of such 
edges, we see that (ii, «2, • • • , *at) has even potential if and only if there exist j±, j 2 , ■ ■ ■ , ]n 
such that (ii, i 2 , . . . ,iff) is sub-positive. 

As one last step, we define the value of the cycle (ii, i 2 , ■ ■ ■ , in), given a choice of 
reaction indices (ji, . . . , j'jv), to be the quotient 
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This expression amounts to an alternating product and quotient of the stoichiometric 
coefficients of the species in the chosen reactions. Its significance will become apparent 
in the next section. 

We are now able to state the main theorem. 

Theorem 3 (Craciun-Feinberg |16j . restated). Suppose an open (homogeneous) CRN has 
the property that if a cycle in G(x) is sub-positive for some chosen reaction indices, then 
it has value 1. If no two sub-positive cycles strongly intersect, then the network cannot 
have multiple equilibria, regardless of the rate constants of the reactions. 

In particular, if no cycles have even potential, then the network is injective. The same 
is true if all stoichiometric coefficients are equal to 1 and no species appears in more than 
two reactions [T6] . 
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The proof can be broken into two components: showing that the cycle condition implies 
a nonzero Jacobian determinant and showing that the nonzero Jacobian determinant 
implies injectivity. 

6.2 The Craciun-Feinberg proof from cycles to the Jacobian 

Theorem 4. Suppose an open (homogeneous) CRN has the property that if a cycle in G(x) 
is sub-positive for some chosen reaction indices, then it has value less than or equal to 1 . 
If no two sub-positive cycles strongly intersect, then the Jacobian det(Jp) is nonzero for 
all x and all positive values of the rate constants. 

Note that we have generalized Craciun and Feinberg's result to include a wider range 
of sub-positive cycles. Naively, we would expect that the new version now covers, in 
addition to the cycles with value 1, roughly half of all other cycles. 

Proof. In the natural expansion of det(Jp), each term is a product of n of the entries J^., 
with one term for each permutaion ii, i 2 , . . . , i n of 1, 2, . . . , n. Craciun and Feinberg, by 
essentially splitting the .1^ into sums of J^, cleverly rearrange the expansion into terms 
corresponding to a permutation z'i, i 2 , . . . , i n and some list ji, j 2 , . . . , j n Each of these 
new terms turns out to be a product of some positive rate constants and concentrations, 
multiplied by det([cj 1 , . . . , CjJ) • det([ej 1 , . . . , ejj), with these matrices of es formed by 
stacking the n row vectors listed@ We would like to show that all such coefficients - the 
products of the determinants - are zero or have sign (—1)". 

Let us consider a nonzero product det([cj 1 , . . . , CjJ) • det([ej 1 , . . . , ejj). In the first 
determinant, there must be a nonzero entry in each row and in each column, so we can 
switch the order of the ih such that the diagonal entries (of both matrices, in fact) are 
all nonzero. Different orderings satisfying this condition correspond to nonzero products 
involving different permutations of the same n reactions ji, j 2 , ■ ■ ■ ,j n - 

Now consider det([ej 1 , . . . , ejj). From our discussion of Theorem [2J we know that the 
nonzero terms in its (natural) expansion correspond bijectively to all sets of disjoint 
cycles in the interaction graph for [ej 1 , . . . , ejj . For this matrix, a cycle consists of 
indices ii,i 2 , . . . ,in such that each ej h i h is nonzero, which is to say that S ih+1 appears 
in reaction j^. Since the diagonal entries are nonzero, we know, as in section IdTT| that Si h 
is on the reactant side of reaction jh for all h, and so Cj h i h and Ji h+1 j h i h are nonzero. 
Thus, the cycles here are exactly the same as those in G(x) such that, not only is Ji h+1 i h 
nonzero, but Ji h+1 j h i h is nonzero as well. Moreover, the signs of the cycles here are the 
same as the original sub-signs, since Cj h i h+1 has the same sign as Ji h+1 j h i h , as observed in 
section E?U Note that if any of these cycles, when regarded in G(x), contain Si h , then 
the outgoing edge from Si h must correspond to reaction jh, meaning that any intersecting 
cycles from [e^, . . . , ej n ] are strongly intersecting in G(x). 

11 Because of the outflow, there are always at least n reactions in the network. 

12 In Craciun-Feinberg, these matrices are given as the transposes of the versions here, but we can make 
the change because taking transposes preserves determinants and flips the directions of all edges, so that 
cycles and their intersections are preserved, with opposite orientation. 
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Suppose we examine a term T in det([ej 1 , ejj) having c = (z'i, 12, ■ ■ ■ , in) 

among its set of corresponding disjoint cycles. This cycle c contributes a factor of 
C = (— l)^" 1 Ylh=i e ihih+i t° the term (where C = T only if N = n), where the 
leading (— l)^" 1 is due to the sign of the permutation associated with c. By the corre- 
spondence of cycles with those of G(x) , we know that if c is sub-negative, then the sign of C 
is (— l) N ~ l (— 1) = (— 1)^, and if c is sub-positive, then the sign of C is (— l)^ -1 . There 
also exists another nonzero term T' in the expansion which differs from T by replacing C 
with C = Ylh=i e jhih> i- e - by taking the identity permutation of the set {it, z 2 , • • • , «tv}. 
Since ej hih = —Cj hih for all h, the sign of C must be (— 1)^. 

If all of the cycles contributing to T are sub-negative, then the sign of T is (— l) n . More 
generally, suppose that for some choice of reactions ji,j2, ■ ■ ■ , jn as above, we have a set Sp 
of sub-positive cycles, all of which are weakly disjoint from each other and have \C'\ > \C\. 
Let Sn be any set of disjoint sub-negative cycles, with total associated product Cn, and 
suppose there are M cycles in Sp that are disjoint from all elements of Sn- If these cycles 
have associated factors Ci, C 2 , . . . , Cm and C[, C' 2 , ■ ■ ■ , C' M as above, then the sum of all 
of the terms in the expansion of the determinant that contain exactly these sub-negative 
cycles and some subset of the M sub-positive cycles is CNYltLiiCk + C' k ), since each 
sub-positive cycle either is or is not present. By the condition \C'\ > |C|, this quantity 
has the same sign as CArllfcliCifc) which is (— l) n . Taking all sets S^, we thus account 
for all of the terms in the expansion of detQe^, . . . , ej n ]), and hence the sign of the 
determinant is (— l) n l^l 

In order to complete the proof of Theorem HI we need to show that det([cj 1 , . . . , Cj n ]) 
is positive. Luckily, the same argument we have just used applies again here, but with 
two modifications. First, the full set of possible cycles is reduced from all of those in G(x) 
to those having all edges sub-negative. Second, the entries in [cj 1 , . . . , CjJ are all positive, 
meaning that the diagonal term in the determinant, and hence the entire determinant, 
has sign +1 instead of (— l) n . 

□ 

6.3 The Craciun-Feinberg proof from the Jacobian to injectivity 

Once we know that all of the terms in the determinant of Jf have the same sign, all we 
need in order for the determinant to be nonvanishing is that some term is nonzero. For 
this, we simply take the product of the diagonal elements, all of which are strictly negative. 
We would then like to prove that the nonzero Jacobian determinant implies injectivity. 
The method here is taken directly from Craciun and Feinberg [T5|, but whereas they only 
deal with CRNs, we will expand the proof to a more general class of functions. 

Let rrii and Vij be nonnegative integers for i — 1, 2, . . . , n and j — 1,2, ... , m,. For real 
numbers A = (aij), let F j4 (x) : (IR + ) n — > M. n be a polynomial function with components 
Ff (x\, X2, • • • , x n ) = ^^=1 a ijX Vii , where we use the shorthand u v = J"J itf . 

13 Craciun and Feinberg [TB] state and prove this result with the requirement that \C'\ = \C\, using 
strict cancellation of terms. 
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From F A , remembering that the Xi are positive, we define the associated Jacobian 
matrix J f a = (Jjfc), where 

dx k ^ x k 

Theorem 5. Fix mi and v;j as above. The following are equivalent. 1: The function F A 
is injective for all A. 2: The determinant of 3 f a is nonzero for all A and all x £ (M + ) n . 

Proof. Suppose that for some values of x and the a^, the determinant of J f a is zero. 
Since J f a is a linear transformation, this means that there is some nonzero vector y £ M n 
such that Jp^y = 0. From above, the vector Jp^y has components 



(J F Ay) l = ( Yl ~~' , '.' xVi ' i 



fe=i \i=i fe 



which we assume are all zero. We claim that these equalities are equivalent to 

mi n 

vijVijkSk = o 

i=i fc=i 

for all z, where rjij = ai-,x Vi j and 5k = Vk/xk- All of these are real numbers, and not all 
of the 5k are zero. In the forward direction, the substitutions are straightforward. In the 
reverse direction, we can obtain the implication by, for example, taking all of the x^ to 
be 1 and then letting = rjij and yk = 5k- 

Next, we claim that this last set of equalities is equivalent to the set 

y~] Kij ■ ^exp ^2 v ijk5k^j - 1^ = 



for all i, where 



Vij ■ Y Vi rth 



k=l 



K = 

exp I Y v ijk5 k ) - 1 

Again, the forward substitution is no problem, with the caveat that the denominators may 
be zero, in which case we just choose arbitrary for those indices, and the equalities 
hold regardless of what we choose. The reverse substitution, for the rjij in terms of the K^, 
is analagous. 

From this set of equalities, we can pass to 
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for all i, taking all the and bk positive and ak/bk = e 5k for each k. In reverse, we simply 
let 8k = ln(afc/fefc) for each k, noting that the condition that not all of the 8k are zero is 
equivalent to the condition that b. Finally, the above are equivalent to the relations 

mi 

^L ir (a v «-b v «) = 

i=i 

for all i, where Ljj = Kij/b Vi i, or equivalently, Kij = L^b Vi j. And, at last, we see that 
our train of equivalences, starting at a determinant of zero for some point in the domain 
and some coefficients, ends with a lack of injectivity for some (different) values of the 
coefficients. This proves the theorem. 

□ 

For the case of a CRN, with reactions of the form Cj\Si + Cj 2 S 2 + ■ • • + Cj n S n — > 
djiS\ + dj 2 S 2 + ■ — h dj n S n , we would have Fj = dxi/dt = Y^jLi e jikj^- Ci , so that rrii = m 
for all i, a.ij = ejikj, and Vij = Cj for each i. Note that in the proof above, we can restrict 
our coefficients to be positive by using the fact that r(e r — 1) > for all real numbers r, 
with equality only when r = 0. 



6.4 Applying the theorem 

As discussed at the beginning of section El the assumption of outflow for CRNs is critical 
to the Craciun-Feinberg theory (and well justified in many cases). Experimental and 
theoretical evidence suggests, though, that closed systems, having at most partial outflow, 
are easier to constrain in their behavior than are open ones [T3], [351 HZ]- Feinberg's 
deficiency theory, for example, gives extensive information for networks with deficiency 0, 
partial information for deficiency 1, and no information for higher values, and closed 
systems tend to have low deficiency [HI EH [22]. Additionally, examples of multistability 
in the laboratory have tended to come from closed systems [15] . 

To illustrate the usefulness of the Craciun-Feinberg theorem, however, let us consider 
the first experimental example of bistability in an open homogeneous system (at least 
according to its publishing author), which was found by Degn [19] in the peroxidase- 
catalyzed oxidation of NADH. He assumes a very simple mechanism, consisting only of 
Michaelis-Menten enzyme action with one additional substrate-inhibition step] 14 ! 

1. E + S^ES 

2. ES^E + P 
3,4. ES + S^ES 2 . 

The form ES 2 is a deactivated enzyme-substrate complex. While we cannot prove that the 
network supports bistability, we can note that, despite its simplicity, the interaction graph 
contains two strongly intersecting sub-positive cycles: E =^ ES =>- E (using reactions 1 

14 The irreversibility of the first step is irrelevant. 
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and 2, both sub-positive edges) and E =^> S =>- ES =>- E (using reactions 1, 3, 2, yielding 
sub- negative, sub-negative, and sub-positive edges). 

By contrast, consider two possible schemes involving product inhibition: 



1,2. E + S^ES 

3. ES^E + P 

4,5. E + P^EP 

or 4', 5'. ES + P ^ ESP. 

Neither of these mechanisms contains any stoichiometric coefficients greater than 1, and 
for neither one does the interaction graph have two strongly intersecting nontrivia^ 
positive cycles. Thus, by Theorem [3j the networks are injective. As a result, we can 
confirm for Degn that the ability of the peroxidase system to generate bistability cannot 
be due to product inhibition, even though all three mechanisms discussed here have very 
similar forms. 



7 Comparing and reconciling the approaches 

The theorems discussed above come from a wide range of fields. Jelonek and Pinchuk 
approach the real Jacobian problem as pure mathematicians, Clarke writes as a chemical 
physicist, Gale and Nikaido are mathematical economists, Thomas and Soule come with 
a biological perspective, and Craciun and Feinberg are chemical engineers. Given the 
similarities in the methods of the Gale-Nikaido/Thomas-Soule and the Craciun-Feinberg 
theorems in particular, though, it is not surprising that the results are closely linked. 

7.1 Introduction to the case of no concentration-ambiguity 

The relationship between the Thomas-Soule and Craciun-Feinberg theorems is particu- 
larly strong when all of the elements in the Jacobian have a fixed sign for all values of the 
state variables (or, in our terminology, no edges in the interaction graph are concentration- 
ambiguous). In general, the assumption of no concentration-ambiguity is a fairly restric- 
tive one, since the entries in the Jacobian can be complicated functions. As we will see, 
though, it has been made in numerous places, albeit without extensive justification. Our 
analysis begins with the following lemma. 

Lemma 6. In the interaction graph of a CRN, an edge is concentration- ambiguous if and 
only if it is reaction- ambiguous. 

Proof. First, suppose we have a CRN such that the edge from S^ to Si 2 is not reaction- 
ambiguous, which is to say that all of the terms Ji 2 ji 1 have the same sign. Since all species 
concentrations are positive, we see immediately that has the same sign for all x, and 
hence the edge is not concentration-ambiguous. 

15 See section [731 
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Next, suppose the edge from to Si 2 is reaction- ambiguous; we would like to show 
that as x varies, is sometimes positive and sometimes negative. Given a nonzero 
term 

/. .. = £z*i e .. /• x «li 

its sign is determined by that of €ji 2 . Looking at reaction j (as in section 16.11) . we see 
that J%2ji\ 

is positive when Si 2 is on the product side and negative when Si 2 is on the 
react ant side. 

Note that the rate term above contains a factor of Xj 2 if and only if Cji 2 > in the 
exponent, i.e. if Si 2 is on the reactant side of reaction j. As a result, in the limit of 
very small values of Xi 2 , the sum of the negative terms in J isil will approach 0, while 
in the limit of very large values, this sum will approach — oo. Meanwhile, the positive 
terms in J i2il are unaffected as Xi 2 varies. By continuity, then, the edge from to Si 2 is 
concentration-ambiguous. 

□ 

This lemma gives us a clear characterization of concentration-ambiguity for CRNs, 
namely that it obtains precisely when a species serves as a direct activator for another 
species Si 2 in one reaction and as an inhibitor of it in another reaction. Based on this con- 
dition, we might indeed expect to find chemical systems without concentration-ambiguity 
reasonably often, although there are certainly well-known cases that are ambiguous, espe- 
cially in gene regulation [TT] . For example, to use we have seen, the first product-inhibition 
system discussed in section 16.41 has no concentration-ambiguity, while the second does. 
The first system also provides another illustration of the importance of assuming that 
enzyme-driven processes are broken down into elementary reaction mechanisms (as in 
section [373]) . since the edge from E to P would be concentration-ambiguous if the first 
two reactions were replaced by the single step E + S — > E + P . 



7.2 Thomas-Soule and Craciun-Feinberg 

Before Soule's completion of Thomas's conjecture, Gouze [26J, Plahte-Mestl-Ohmolt |42j . 
and Snoussi [IE] all published calculus-based proofs of the conjecture with the additional 
assumption of no concentration-ambig uity@ The Craciun-Feinberg theorem also gives us 
a simple proof of the conjecture in this follows. Suppose that for a given CRN, 

G(x) has no positive cycles and no concentration-ambiguous edges for any x. The sub-sign 
of any cycle will always be the same as its sign, i.e. it will be negative. Thus, Theorem [3] 
immediately tells us that, as long as our system is open, it must be injective. In fact, 
we obtain the stronger result that injectivity holds regardless of the rate constants. Our 
conclusion is most similar to that of Snoussi, since he also assumes (in effect) an outflow 
or degradation term for each state variable [4~8] . 

We can also make some (but not complete) progress with this argument when edges 
are allowed to be concentration-ambiguous. Again, we would like to show that negative 

16 The result of Cinquin-Demongeot is somewhat more general. 
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cycles must be sub-negative, so that we can apply Theorem[3j Suppose C = (z'i, 12, ■ ■ ■ , ijv) 
is a negative cycle that is sub-positive for some choice of reactions ji,j2, ■ ■ ■ ,3n- There 
must be at least one edge in C that is reaction-ambiguous, or else C would be positive as 
well as sub-positive. Without loss of generality, let some terms with 

positive signs and some with negative signs. 

Let us vary Xi 2 , as in the proof of Lemma [6], until there is a change in the sign of J^n- 
Unless the sign of another edge in C changes simultaneously, we will have a contradiction, 
because C will be positive for slightly larger or smaller Xi 2 . In other words, at any 
concentration vector such that = 0, we must have Ji h+1 i h = for some other edge in 
the cycle. This gives us a finite set of polynomials such that all of the zeroes of one of 
them (J^ = 0) that are in (R + ) n are zeroes of at least one of the others. This situation 
is mathematically possible, but it seems physically unrealistic, a "measure-zero" case. 

It is true, however, that this case could arise. For an example, we turn back to the 
CRN in section 15.21 assuming outflow reactions for each species in addition to the six 
reactions shown. As noted above, the cycle C^A=>Z=>D=>B=>C is always 
negative, yet it is sub-positive for the correct choice of reaction indices. This example is 
also a good illustration, however, of our statement that this is a "degenerate" or measure- 
zero condition: mathematically, there may be nothing exceptional about the system, but 
in real life, we could not expect any pairs of rate constants to be precisely equal. 

7.3 Monotone systems and Polya's permanent problem 

It is worth taking a moment to acknowledge two other sources of related results. The first 
is the theory of monotone systems, as discussed in the work of Angeli and Sontag (2j[3j,|4]. 
A CRN gives rise to a monotone system provided that none of the edges in its interaction 
graph are concentration-ambiguous and all of the cycles are positive. One would expect 
these conditions to rule out a rather large class of networks, but when they apply, the 
theory gives more extensive information about sufficient conditions for multistability and 
the locations of steady states than do any of our other approaches. The example of 
MAPK cascades figures prominently in their papers; see section [HJ for a discussion. A 
special subcategory of monotone systems is that of cooperative systems, which are those 
with interaction graphs having all edges positive for all x [2] . 

Polya's permanent problem asks when it is possible to relate the determinant of a 
matrix to its permanent, which is formed by adding the terms in the standard expansion 
of the determinant without their associated permutation-signs attached. This question 
is not directly one of injectivity, but as in section 16. 2\ it is very much related to the 
interplay between matrices with signed elements and their related graphs (6j. In fact, the 
problem has appeared in a number of contexts and forms, even including the calculation of 
resonance structures in organic chemistry [32J. Combinatorial extensions of the problem 
have been treated extensively in the graph theory literature [3TI I3"8] , especially with regard 
to conditions for the existence of various cycles in a number of different kinds of graphs. 
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7.4 Craciun-Feinberg and Gale-Nikaido 



We have already seen how Craciun and Feinberg's theorem can imply that of Thomas- 
Soule in the case of CRNs. Here, we show how Gale-Nikaido can be used to prove the 
second half of Craciun-Feinberg. 

Suppose that the cycle conditions of Theorem[3]are met for a given CRN, implying that 
all of the determinants det([ej 1 , . . . , ejj) have the same sign. This fact depends only on 
the network's structure and not at all on its rate constants. Now, consider any subgraph 
of the interaction graph corresponding to a subset of N of the species. Any cycles in 
this subgraph are also present in the full graph, which means that the hypotheses of 
the theorem are still satisfied, and hence there is no way to generate coefficients with 
a sign of (— 1) N ~ X in the principal minor corresponding to this subset. (Because of the 
diagonal term, the minor is strictly nonzero.) Thus, the matrix J_f is a P-matrix, and 
by Gale-Nikaido, the full network is injective. 

This proof does not use Gale-Nikaido to its full potential, since the hypotheses of 
Gale-Nikaido do not require that minors be nonzero for all values of the rate constants. 
Ideally, one would be able to find conditions, cycle-based or otherwise, that would give 
nonvanishing minors for a specific set or region of rate constants, even for a network that 
has some "bad" cycles in it under Craciun-Feinberg. 

7.5 A note on reversible reactions and related issues 

While the Thomas-Soule theorem gives an elegant necessary condition for multistability, 
its power in the case of CRNs is limited by the prevalence of positive cycles. Any reversible 
reaction, for example, even A ±=p B, gives rise to a positive cycle, with positive edges from 
A to B and vice versa. Regardless of the outflow properties, though, this reaction only has 
one equilibrium. Even a single reaction A+B — > C, for that matter, yields a positive cycle, 
this time with negative edges between A and B. These "trivial" examples might make us 
pessimistic about the applicability of any of our theorems, but, luckily, Craciun-Feinberg 
provides some relief. 

Let us assume degradation or outflow of all species, so that the Craciun-Feinberg 
results apply. We have seen in section [6^2] that we can rearrange the expansion of det(Jp) 
such that the coefficient of each term is of the form det([cj 1 , . . . , Cj n ]) • det([ej 1 , . . . , ejj). 
As long as the signs of these coefficients are all the same, the CRN will be injective. If we 
form a cycle with any repetition among the indices ji,j2, ■ ■ ■ ,3Ni then the first of these 
determinants will be zero, since two rows will be identically Thus we can declare that 
the A + B — > C issue in the previous paragraph can be ignored. 

With regard to reversible reactions, Craciun and Feinberg combine both directions into 
a single reaction-node in their graphs, which eliminates the possibility of using a reaction 
and its reverse within the same cycle. Within our framework, we can justify simply 
ignoring such cycles by noting that they lead to the presence of two rows in the second 

1 In the context of their theorem, this issue never arises, due to the strict assignments of reaction 
indices jh (see section \672\ . 
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matrix above that differ by a factor of -1 throughout (and are thus linearly dependent). 
Graphically, we can incorporate this exclusion by drawing single, two-sided arrows for 
reversible reactions in the interaction graph. 

Of course, these trivial examples are not the only ways in which one of the determinant- 
coefficients above may vanish due to linear dependence. All that is necessary is the 
alignment of a few stoichiometric coefficients. In general, though, this will be hard to 
detect without checking individual cycles by hand. 

As a final note, while Craciun and Feinberg's statment of their theorem eliminates 
these two forms of cycles, it simultaneously introduces at least two other potential trivi- 
alities. First, their main theorem uses the non-oriented SR graph, some of whose cycles 
might not be true, oriented cycles. Second, their method for drawing edges in the SR 
graph allows false edges between species that both appear on the product side of a given 
reaction, for example in the "cycle" A =>- B =>- A arising from the reaction C — > A + B. 

8 Multisite protein phosphorylation 

Perhaps the most-studied example of a bistable system in biochemistry is the eukaryotic 
cell-signaling module based on mitogen-activated protein kinase (MAPK) j2j [3J HI 0, [23J 
1271 1281 [36] . A kinase is an enzyme that acts by attaching a phosphate group, usually 
taken from a molecule of ATP, to another molecule in order to effect a change in that 
molecule's chemical behavior. An enzyme that removes a phosphate group is referred 
to as a phosphatase. When abstracted, the key feature of the MAPK system is that 
MAPK is itself regulated by phosphorylation, and moreover, that it has two distinct 
phosphorylation loci. Let us write Si for an z-times phosphorylated molecule of MAPK, 
for i = 0, 1, 2; E and F for the kinase and phosphatase that act on the Si, and ESi 
and FSi for the various intermediate enzyme-substrate complexes. Part of the beauty of 
the MAPK system is that it can be accurately modeled using nothing more than repeated 
applications of Michaelis-Menten enzyme action [7J [27J [28], [36] • Thus, the full mechanism 
is the following. 
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Many of the authors cited above have noted the connection between MAPK cascades 
and positive feedback, and indeed, we can see that the network contains a multitude of 
positive and sub-positive cycles, at all levels of the structure. Each of the four phospho- 
rylation and dephosphorylation events (reactions 1-3, 4-6, 7-9, and 10-12) contains its 



23 



own two-element sub-positive cycle, for example E =>■ ES$ E using reactions 1 and 3. 
We can also split the process into events at the first site (reactions 1-3 and 10-12) and 
those at the second site (4-9), and here we find the cycles So =>• ESo =>• Si =>■ FSi =>- So 
(reactions 1, 3, 10, 12) and Si =>- ES± =>- S 2 =>• FS^ =>• Si (reactions 4, 6, 7, 9). 
These last two are significant because they strongly intersect with the sub-positive cycles 
S Q => E => S!=> F=> S (reactions 1, 4, 10, 12) and S 2 => F => Si => ESi => S 2 (re- 
actions 7, 10, 4, 6), respectively. While we cannot use Theorem [3] to prove multistability 
based on these intersecting cycles, we can see, given that multistability is present, that 
they are responsible. 

By contrast, consider a protein that is phosphorylated at only one site, with a mech- 
anism consisting of reactions 1-3 and 10-12. There are still positive cycles, so we cannot 
apply Thomas-Soule, but now, the last two cycles listed above no longer exist, and in 
fact, there are no longer any strongly intersecting sub-positive cycles. Thus, by Craciun- 
Feinberg, the one-site system is injective. 

9 Concluding remarks 

Multistability theory can be both delicate and forgiving. Through an investigation of 
several approaches to the theory, we have seen that, just as multistable systems can be 
perturbed from one equilibrium state to another, slight differences in structure can push 
a network from the realm of the injective to the realm of the multistable. At the same 
time, it is remarkable that the theorems we have examined, in particular that of Craciun 
and Feinberg, allow us to conclude so much from so little. 
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